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Abstract- The various scientific laws, principles are used to develop 
mathematical models. Often the differential equations are used to describe the 
system. But, formulating the differential equations for most problems is difficult 
and hence obtaining the solutions by exact methods of analysis is a formidable 
task. The present research discusses the issue of the finding the field variable 
deviation for quadratic area of one dimensional continuum. The analysis 
provides the percentage error in the field variable for the selected trial functions. 
Approximate method is useful, if it is integrated with computer for the problem 
involving a number of complexities without making drastic assumption which 
otherwise complicated to attempt by classical methods. A genuine necessity for 
obtaining precise solution for the different numerical approximation methods is 
overcome by developing in-house computer program. Graphically results can be 
displayed to know the effect of considered weights and the constants assumed. 
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I. INTRODUCTION 


The model problem considered is of an axially 
loaded bar having quadratic function of area. The unknown 
variable is the axial displacement of one dimensional 
continuum, u(x) is attempted in the present research by means 
of numerical analysis technique where the basic inputs to a 
problem are known with arbitrary basic data. The models are 
often described in terms of algebraic, differential and integral 
equation which relates quantities of interest. Mathematical 
model is a set of relationship among variables that are 
expressed in essential features of physical system or the 
process in analytical terms. As a result, physical system or 
process can be described by means of mathematical model. 
The relationship that governs the system may take the form 
of algebraic, differential and integral equations. 

The governing equations and boundary conditions 
are useful means to express the mathematical model. In 
general, engineering problems is often needs to be addressed 
considering the domain of interest. The differential equations 
are employed to illustrate the physical significance of the 
system under the domain of consideration. For such system, 
analytical solution can be found out using differential 
equation and boundary conditions. To establish mathematical 
model it is required to describe the relationships of primary 
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and secondary variables. The solution characteristics can be 
expressed in the algebraic form. For the well defined simple 
engineering problems, there are standard known analytical 
techniques available to find solution. But, there are many 
engineering problems for which exact solution cannot be 
easily obtained using analytical techniques. The failure to 
obtain an exact solution may be attributed to either the 
complex governing differential equation or the difficulties in 
dealing with the boundary and initial conditions. To deal with 
such problems, numerical approximation is best alternative. 
An approximate technique may become the effective tool to 
define the nature of unknown variable and enable to 
determine the solution at discrete points in the problem 
domain. 

In the present paper the numerical techniques like 
Direct Integration, Method of Weighted Residuals, Ritz 
Method and FEM are used to determine the deviation in the 
field variable at the selected node position. For this purpose 
the one dimensional continuum is selected as shown in Fig. 1. 
The computer programming technique is employed to 
minimize computational cost and time. The developed 
computer program uses MATLAB tool and developed sets of 
subroutines and functions. The programs are intended as a 
primary component for generalization to obtain the results at 
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any position on the continuum. The applications can be 
extended for variety of similar complex problems. 


Il. DEVELOPMENT OF GOVERNING EQUATIONS 


Consider an axially loaded 1D bar having varying area A 
with one end fixed at x=0 and externally applied force F at 
x=L as shown in Fig.1.The free body diagram (FBD) is 
shown in Fig. 2. The axially loaded continuum considered 
here has length L=2 m and area as exponential function of x 


given by A=A,e ©, E=200 GPa and F=1KN. 


F= 10008 


E=200 GPa 
AgeSé10" wr 





Fig.1 Axially loaded 1D continuum 
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Fig.2 Free body diagram 


Governing equation is developed for 1D continuum having 
varying area. The unknown variable is the axial displacement 
given by u(x). 

From the free body diagram (FBD) and using equilibrium 
equation, 

oA += (cA)dx — cA =0 (1) 
Where, A indicates the area which is assumed to be a 
continuous function with respect to x. The equation (1) 
becomes, 


<- (cA) = 0 (2) 

From Hooke’s Law, within elastic limit, Stress « Strain, 
du 

orao=E e (3) 


Where E is the modulus of elasticity of the system, 
substituting o in equation (2) 
The governing equation is, 
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< (EAS) =0 , 0<x<L (4) 
There are two boundary conditions associated with the 
problem, 

The essential boundary conditions is. 


u(x = 0) =0 (5) 
and the natural boundary condition, 

ou a tee 

gt = Se (6) 


Where F is the axially applied load 
I. DEVELOPMENT OF ROUTINE SCRIPT 


MATLAB provides interactive platform for executing 
numerical computations. The program code is written in 
MATLAB 7.10.0 (R 2010a) version. In-house computer 
program is developed different approximation numerical 
methods to determine the solutions to overcome long and 
tedious task which otherwise need to be attempted by hand 
calculation. The important steps which are used for writing 
the program script are indicated in the block diagram shown 
in Fig3. 


Execution of 
Mathematic 
al operations 
according to 
selected 
method 


Selection of Selection 
problem and 
indentifying 
corresponding 
boundary 
conditions 


Developm 
ent of of trial 


Governing function 


equation 


Formulation 


Plot Finding and 


Displace Displace 
ment ment 


Finding 
constants 


execution of 


simultaneou 


of trial s equations 
function function function 


results 





Fig 3. Block diagram representation used for developing 
program. 


IV. DIRECT INTEGRATION APPROACH 


The direct integration technique is analytical way of finding 
exact solution [8]. The result of approximate method is 
verified using classical methods i.e. Direct Integration 
method. The technique is used for obtaining the result and 
finding the error of field variable at the discrete points. 
Applying boundary condition, u(x = 0) = 0 

Substituting values of F, E and A at x=2 and applying 2" 


boundary condition 
F 


Equation 6 becomes 8 (7) 
dx AE 
Integrating the governing equation (4),(EA “) =c(,> 
oso (8) 


dx EA 
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du = —*—, dx (9) 


ExAoxe T 


* 
1 C1xel 


Again Integrating, u = eee + C2 (10) 


Now, applying 2" boundary condition to equation (8) i.e 
substitute x=2, the result is C, = 1000 
Put the value of C, in equation (10) at x=0, 

(.==-2% 10" 
Substituting the value of c, and cz in equation (10), the 
generalized equation for displacement is, 


2x1000*100 _s5 
200x109x5x10-4) — 2x 10 (11) 


The displacement at interval of 0.25 m increment position 
from the starting point is evaluated and depicted in Table 1 
and the graphically results are represented in Fig 4. 





u(x) = 


A. Program code for Direct Integration method 

The developed program calculate‘cl’ and ‘c2’, the constants 
of the governing differential equation. The deformation is 
calculated in three statements of program. ‘X’ is an array 
having different values of distances from fixed end. The 
second step call the value of X and final deformation values 
are obtained in the third step. 


F=input('Enter load acting on body in "N":-'); 
E=input('Enter the value of Youngs modulus in "GPa":-'); 
E=E*(10‘9); % to convert in GPa 
A=(0.0005)*(exp(-x/2));% Cross-sectional Area 
x1=input('The starting position of deformation in '"m":-'); 
x2=input('The ending position of deformation in '"m":-'); 
n=input('The number of parts the body is to be divided:-'); 
X=x1:(x2/n):x2; % produces array of distance from fixed end 
a=5*(10%-4); 

D=(R/(E*(a*(exp(-1)))));% First differential of u wrt x 
cl=(D*(E*(a*(exp(-1)))));% Constant of integration | 
c2=-(4*cl *(exp(x1/x2))/(E*(a)*x2));% Constant of 
integration 2 

k1=(4*(exp(X/x2))*c1)/((E)*x2*5*(104-4)); 
ue=k1+c2;% Equation of deformation 


V. APPROXIMATION METHODS 
Approximate methods use numerical technique to obtain the 
solution from the governing equations. The present section 
discusses the issue of the finding the field variable deviation 
for the various positions on the continuum. The different 
standard numerical techniques are applied to these one 
dimensional continuum systems to determine results at 
different nodes. 


A, METHOD OF WEIGHTED RESIDUALS 

The weighted residual methods are based on the assumption 
of an approximate solution for the governing differential 
equation. The assumed solution must satisfy the initial and 
boundary conditions of the given problem. As the assumed 
solution is not exact, substitution of the trial solution into the 
differential equation will lead to some residuals or errors. 
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Each residual method requires the error to vanish over some 
selected intervals or at desired points called nodes. 

The governing equation for the above considered problem is 
equal to residual(R) when wu = Uapproximate, Where 
Uapproximate iS a polynomial function which describes the 
solution at the discrete points [6]. 

The governing equation for these methods is 


d duapprox\ _ 
< (EA Severe) = R (12) 
Considering trial approximate solution, 
U = Uapproximate Ao + A,X + px? (13) 
* ppt =a,+2a,x (14) 
Applying first B.C. u=0 at x=0 
“ ao = 0 
Therefore the equation (13) becomes, 
U = UapproximateA1X + px? (15) 


This function is subsequently used for further calculation. 
The weighted residual methods uses virtual work principle 
and the integrals error function are defined as 
ai, WR a0, 1128 ca. (16) 
L d d 

+ fy Wi x £ (EAS) dx = 0 
The equation is integrated by parts 

du du L 
# ((WAE = Jyar — (WiAES) x20} — Jy EA 
(17) 
Where W; are the weights, which are different for each 
method. 
The constants of trial solution a, and az are calculated by 
solving equation (12) for different weights. The deformation 
is calculated by substituting them in equation(15) 

“. U = Uapproximate 21% + nx? 


du dW; 
dx dx 


dx =0 


In the developed program the constants are’al’and ‘a2’ 
are declared. The letter ‘x’ is variable in ‘u’ of trial solution 
assumed. The first order differentiation of ‘u’ is calculated as 
‘du’. The weights are defined as ‘wl’ and’w2’ which will 
modify as per the methods considered for analysis. The 
differentiation of these weights is calculated as ‘dwl’ and 
‘dw2’. The differential equation ‘fl2’ is formulated by 
building the equation in terms of ‘f11’. Similarly the second 
equation ‘f22’ is formulated and integrated over the domain. 
The integration provides us with two simultaneous equations 
‘cl’ and ‘c2’. These equations are solved simultaneously and 
their values are saved in ‘a’, which is an array holding the 
solution. The constants of trial solution are retrieved from it 
and the trial solution is calculated for different values of ‘X’ , 
the distance of a point from fixed end. The common syntax 
applicable for the different Method of Weighted Residuals is 
depicted below with the appropriate comments with each 
statement. 
syms x; 

A=(0.0005)*(exp(-x/2)); % Cross-sectional Area 
syms al a2; 
u=(al *x)+(a2*(x‘2));% Trial solution 
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du=diff(u,x);% First differential of u wrt x 

wl= %**% Weight defined by particular method 
dw1=diff(w1,x);% First differential of w1 wrt x 

fl 1=diff((w1*F),x);% Part 1 of differential eq(1) 
f12=(dw1*A*E*du);% Part 2 of differential eq(1) 

cl=int(fl 1,x,x1,x2)-int(f12,x,x1,x2); % Total differential 
eq(1) 

w2= %**% Weight defined by particular method 
dw2=diff(w2,x);% First differential of w2 wrt x 
f21=diff((w2*F),x);% Part 1 of differential eq(2) 
f22=(dw2*A*E*du);% Part 2 of differential eq(2) 

c2=int(f2 1,x,x1,x2)-int(f22,x,x1,x2); % Total differential 
eq(2) 

a=solve(cl,c2,al,a2); % Solving eq (1) & (2) simultaneously 
aal=a.al; % Constant | of Trial solution 

aa2=a.a2; % Constant 2 of Trial solution 

uu=(aal *X)+(aa2*(X.%2));% Equation of Trial solution 


The above syntax will remain same only with the changes in 
the weight used corresponding to the method of selection 
considered for execution. 


a) Galerkin Method 

Galerkin Method is a Weighted Residual Methods which 
derives its weight functions from the trial solution. 

Weight function is given by 


W, = “approximate f= 12,3... (18) 
Therefore W, = ei ppreeaate & 
da, 
W, = dUapproximate 
2 daz 
Integral of the error function is Ie Wi.R.dx =0 [2]. 


Considering the weight functions as W,& W,, and finding 
the constants, the trial solution is obtained. 

The weights for this method can be calculated in form of 
wl_gal=diff(u_gal,al_gal);% Weight(1)obtained from u 
w2_gal=diff(u_gal,a2_gal);% Weight(2)obtained from u 


The constants are determined and field equation obtained is, 
Uapproximate = 8.9158 x 10~°x + 4.0871 x 10~°x? 
Displacements at various points are determined and 
represented in Table 1 and the graphically results are 
represented in Fig 4. 


b) Petrov-Galerkin Method 

Petrov-Galerkin Method works on similar approach 
which uses different pairs of weight functions such as 
W, = 1&x, x&x?,1&x? ete. 
Considering the weight functions as x*& x? for further 
calculation. 
The weights for this method can be written in program as 
wl_petro_gal=x%2; % First Weight considered 
w2_petro_gal=x%3; % Second Weight considered 
Again the same integral error function is used which is 
symbolized as 
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Integral of the error function is de Wi. R.dx =0 [2]. 
The constants are determined and field equation obtained is, 
u = 8.9158 x 10°-°x + 4.0871 x 10°°x? 

Displacements at various points is depicted in Table 1 and the 
graphically results are represented in Fig 4. 


c) Sub domain Method 
The domain of Fig.1 is divided in two halves and the 


analysis is carried out with single weight function [7]. 

a : dw 
Considering weight as w=x,  *. ae rl 
The domain is 0 to 2m which is divided into two halves from 
0 mto 1 mand | mto 2 m. The governing equation (17) used 
in this method which is multiplied by considered weight 


function. The equation is expressed as 
J, (W x R)dx = 0 (19) 
For sub domain 0 to = and 7 to L equation (19) gives two 


equations which when solved simultaneously provides the 
constants of trial solution a, and a2. Additional variable in 
term of x3=(x2-x1)/2; is defined to indicate intermediate 
limit of integration. 

To execute the program following changes are carried out. 
cl=int(fl 1,x,x1,x2)-int(f12,x,x1,x3); % eq(1) 

c2=int(f2 1,x,x1,x2)-int(f22,x,x3,x2); % eq(2) 

The constants are determined and field equation obtained is, 
Uapproximate—9 X 10~°x + 4.09 x 10°°x? 

Displacements at various points is depicted in Table 1 and the 
graphically results are represented in Fig 4. 


d) Least square Method 

The Least Square Method is a weighted residual method 
which derives its weight functions from the Residual itself 
[18]. Weight function is given by 


dR ne 
W; = Ga: i= 1,2,3... (20) 
dR dR 
Therefore W, = ae & W= ae, 


Weights are multiplied with Residual and integrated over the 
domain [2]. 

The weights for this method can be calculated in form of 
res=diff((E*A*du_les),x); % Residual 
wl_les=diff((res),al_les); % Weight (1) obtained from 
Residual 

w2_les=diff((res),a2_les); % Weight (1) obtained from 
Residual 

The constants are determined and field equation obtained is, 
Wapprodmate = 90-1946 X10" °x +3.9497 x 10° 
Displacements at various points is depicted in Table 1 and the 
graphically results are represented in Fig 4. 


B. RITZ METHOD 

In continuum mechanics, a system can be described in 
terms of an "energy functional", which measures the energy 
of proposed configuration. It is often impossible to analyze 
all of the infinite configurations of system with the least 
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amount of energy; an approximate numerical computation is 
essential [7]. Ritz method is applied to the above discussed 
static equilibrium problem. The trial functions selected above 
is substituted in to the functional. The 72 is equal to 
summation of strain energy and applied work potential. 
Functional governing equation is 


Li du\ 
m= fo >x AE (S) dx — Uyer_y X F} (21) 


Now differentiate equation (21) with constants of trial 
function a, and az, 


On On 
ia; = 9 (22) and = =0 (23) 


Solving simultaneous equation 22 and 23 the constants are 
determined and field equation obtained is , 

Uapproximate = 1.534 x 107°x + 1.111 x 107®x? 
Displacements at various points is depicted in Table 1 and the 
graphically results are represented in Fig 4. 


a) Program code for Ritz method 

The program calculates the deformation using Ritz method. 
The approaches are similar to MWR. The variables which are 
required to define formulated for the trial solution are 
selected. The 7 function is calculated by evaluating ‘p1’, ‘p2’ 
and ‘ulr’. Two equations are formed by differentiating the 
‘pi’ function with respect to the constants of trial solution. 
The simultaneous equation is solved to calculate the value of 
‘al’ and ‘a2’. 

syms x al a2 

u=(al *x)+(a2*(x‘2));% Trial solution 

du=diff(u,x);% First differential of u wrt x 
pl=(0.5*E*A*(du‘2));% Part 1 of Pi inside function 
ulr=((al *x2)+(a2*(x2%2)))*R;% Part 2 of Pi function 
p2=int(p1,x,x1,x2);% Part 1 of Pi final function 
pi=p2-ulr;% Total Pi function 

d1pi=diff(pi,al);% First differential of pi wrt al eq(1) 
d2pi=diff(pi,a2);% First differential of pi wrt a2 eq(2) 
a=solve(d1pi,d2pi,al,a2);% Solving eq(1) & 
simulteneosly 


eq(2) 
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Fig 4 Displacement plot for the considered continuum 


Above table 1 indicates the comparison of results 
between different method. Out of which we can conclude 
“Lest Square Method’’ to be more accurate on observing the 
above table we come to know that the result of first two field 
variable are more precise while there is more deviation of the 
result in the intermediate phase, mean while the result of 
“Galerkin Method’’ are accurate but for few field variables 
and again for the extreme phase we get precise result for 
“Lest Square Method’’. 

The percentage error is calculated with respect to direct 
integration method and it is represented in Table 2. 


Table 2 Percentage error in deformations 






































; . Methods | Galerkin | Petrov- Sub- Least Ritz 
al=a.al;% Constant | of Trial solution —> Galerki | domain | square 
a2=a.a2;% Constant 2 of Trial solution : n 
uu=(al *X)+(a2*(X.%2));% Equation of Trial solution Distance 

% : : ‘ : s(m) 0 0 0 0 0 
Executing the program the displacements at various points is Om 
depicted in Table | and the graphically results are represented 0.25 0.1541 | 0.146 | 0.0639 | 0.041 | 0.154 
in Fig 4. 62 1 35 1 
Table 1 deformation results for the numerical methods 0.5 0.0299 | 0.100 | 0.0334 | 0.021 | 0.029 
considered 35 5 1 9 
Methods | Direct Galerkin | Petrov- Sub- Least Ritz 0.75 0.0055 0.064 0.0088 -0.003 0.005 
istance | integration Galerkin | domain | square 91 5 
(am) sige snes. || scaaee sane | ie 1 -0.0077 | 0.038 | -0.0061 | -0.013 | -0.007 
x 10-6 1.25 - 0.021 | -0.0138 | -0.017 | -0.014 
0.01498 3 
0 0 0 0 0 0 0 1.5 -0.0148 | 0.011 | -0.0147 | -0.015 | -0.014 
0.25 2.66 2.25 2.27 2.49 2.55 2:25 
iG 5 68 55] S11 5.49 556 5 5] 1.75 -0.0122 | 0.005 | -0.0125 | -0.010 | -0.012 
0.75 9.09 9.04 8.50 9.01 9.12 9.04 2 0.00146 | 0.012 0.0008 0.005 0.001 
1 12.97 13.07 12.47 13.05 | 13.14 | 13.07 
1.25 17.36 17.62 16.99 17.60 | 17.66 | 17.62 
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VI. CONCLUSIONS 

The issue of selecting the numerical methods to obtain 
the approximate solution for the system is addressed with 
different numerical approximation methods. The weighted 
residual methods explicitly include the errors in the field 
variable which is found out using the developed program. 
Fine variation in the field variable is possible with the 
developed program by using the desired steps of increments. 
The best convergence is found out by executing these 
numerical techniques to select the appropriate regions of 
operations. The smooth convergence close to the exact 
solution is essential for a best possible selection of trial 
solution which is possible with the developed program. 
Selection of trial function can be forecasted based on the 
results of field variable error percentage. The computer 
programming technique is employed to minimize 
computational cost and time. 
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